Molecular Dynamics Simulation on a Metallic Glass Ni .2Zr .8-System: Non-Ergodicity 
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At the present paper we have computed non-ergodicity paramater from Molecular Dynamics 
(MD) Simulation data after the mode-coupling theory (MCT) for a glass transition. MCT of dense 
liquids marks the dynamic glass-transition through a critical temperature T c that is reflected in the 
temperature-dependence of various physical quantities. Here, molecular dynamics simulations data 
of a model adapted to Nio.2Zro.8 are analyzed to deduce T c from the temperature-dependence of 
corresponding quantities and to check the consistency of the statements. Analyzed is the diffusion 
coefficients. The resulting values agree well with the critical temperature of the non-vanisihing 
non-ergodicity parameter determined from the structure factors in the asymptotic solution of the 
mode-coupling theory with memory-kernels in "One-Loop" approximation. 
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I. INTRODUCTION 

The transition from a liquid to an amorphous solid that 
sometimes occurs upon cooling remains one of the largely 
unresolved problems of statistical physics 0, • At the 
experimental level, the so-called glass transition is gener- 
ally associated with a sharp increase in the characteristic 
relaxation times of the system, and a concomitant depar- 
ture of laboratory measurements from equilibrium. At 
the theoretical level, it has been proposed that the tran- 
sition from a liquid to a glassy state is triggered by an 
underlying thermodynamic (equilibrium) transition Q; 
in that view, an "ideal" glass transition is believed to oc- 
cur at the so-called Kauzmann temperature, Tk- At Tjc , 
it is proposed that only one minimum-energy basin of 
attraction is accessible to the system. One of the first ar- 
guments of this type is due to Gibbs and diMarzio , but 
more recent studies using replica methods have yielded 
evidence in support of such a transition in Lennard-Jones 
glass formers 0, IM 0- These observations have been 
called into question by experimental data and recent re- 
sults of simulations of polydisperse hard-core disks, which 
have failed to detect any evidence of a thermodynamic 
transition up to extremely high packing fractions Q • One 
of the questions that arises is therefore whether the dis- 
crepancies between the reported simulated behavior of 
hard-disk and soft-sphere systems is due to fundamental 
differences in the models, or whether they are a conse- 
quence of inappropriate sampling at low temperatures 
and high densities. 

Different, alternative theoretical considerations have 
attempted to establish a connection between glass tran- 
sition phenomena and the rapid increase in relaxation 
times that arises in the vicinity of a theoretical crit- 
ical temperature (the so-called "mode-coupling" tem- 
perature, Tmct), thereby giving rise to a "kinetic" or 
"dynamic" transition |8|]. In recent years, both view- 
points have received some support from molecular simu- 



lations. Many of these simulations have been conducted 
in the context of models introduced by Stillinger and 
Weber and by Kob and Andersen [9j; such models 
have been employed in a number of studies that have 
helped sha pe our current views about the glass transi- 
tion UHlHEmH. 

In the full MCT, the remainders of the transition and 
the value of T c have to be evaluated, e.g., from the ap- 
proach of the undercooled melt towards the idealized ar- 
rested state, either by analyzing the time and tempera- 
ture dependence in the ^regime of the structural fluctu- 
ation dynamics |l5l Il6l Il7| or by evaluating the temper- 
ature dependence of the so-called g m -parameter 0, 0] . 
There are further posibilitics to estimates T c , e.g., from 
the temperature dependence of the diffusion coefficients 
or the relaxation time of the final a-decay in the melt, as 
these quantities for T > T c display a critical behaviour 
\T — T"c| ±7 - However, only crude estimates of T c can be 
obtained from these quantities, since near T c the critical 
behaviour is masked by the effects of transversale cur- 
rents and thermally activated matter transport, as men- 
tioned above. 

On the other hand, as emphasized and applied in 
pfl l2ll l22| . the value of T c predicted by the idealized 
MCT can be calculated once the partial structure fac- 
tors of the system and their temperature dependence are 
sufficiently well known. Besides temperature and particle 
concentration, the partial structure factors are the only 
significant quantities which enter the equations of the 
so-called non-ergodicity parameters of the system. The 
latter vanish identically for temperatures above T c and 
their calculation thus allows a rather precise determina- 
tion of the critical temperature predicted by the idealized 
theory. 

At this stage it is tempting to consider how well the 
estimates of T c from different approaches fit together and 
whether the T c estimate from the non-ergodicity parame- 
ters of the idealized MCT compares to the values from the 
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full MCT. Regarding this, we here investigate a molec- 
ular dynamics (MD) simulation model adapted to the 
glass-forming Nio.2Zro.8 transition metal system. The 
Ni^Zri-^-system is well studied by experiments p3l l24j 
and by MD-simulations EE EH ES El E3> as it is a 
rather interesting system whose components are impor- 
tant constituents of a number of multi-component 'mas- 
sive' metallic glasses. In the present contribution we con- 
sider, in particular, the x = 0.2 compositions and con- 
centrate on the determination of T c from evaluating and 
analyzing the non-ergodicity parameter, and the diffusion 
coefficients. 

In the literature, similar comparison of T c estimates 
already exist EH EJ E2 for two systems. The stud- 
ies come, however, to rather different conclusions. From 
MD-simulations for a soft spheres model, Barrat et.al.pfT] 
find an agreement between the different T c estimates 
within about 15%. On the other hand, for a binary 
Lennard- Jones system, Nauroth and Kob E2 get from 
their MD simulations a significant deviation between the 
T c estimates by about a factor of 2. Regarding this, 
the present investigation is aimed at clarifying the sit- 
uation for at least one of the important metallic glass 
systems. Our paper is organized as follows: In ScctionlTTl 
we present the model and give some details of the com- 
putations. Section IIIII gives a brief discussion of some 
aspects of the mode coupling theory as used here. Re- 
sults of our MD-simulations and their analysis are then 
presented and discussed in Section llVI 



III. THEORY 

In this section we provide some basic formulae that 
permit calculation of T c and the non-ergodicity parame- 
ters fij(q) for our system. A more detailed presentation 
may be found in Refs. EH El El El EH The central 
object of the MCT are the partial intermediate scatter- 
ing functions which are defined for a binary system by 
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where 



q ■ r 



2=1,2 



(1) 



(2) 



is a Fourier component of the microscopic density of 
species i. 

The diagonal terms a = (3 are denoted as the incoher- 
ent intermediate scattering function 



1 ffi 

F*(q,t) = F £ <exp(zq ■ [r* Q (t) - i£(0)])) 



(3) 



a=l 



II. SIMULATIONS 

The present simulations are carried out as state-of-the- 
art isothermal- isobaric (N,T,p) calculations. The New- 
tonian equations of N — 648 atoms (130 Ni and 518 
Zr) are numerically integrated by a fifth order predictor- 
corrector algorithm with time step At — 2.5xl0 _15 s in 
a cubic volume with periodic boundary conditions and 
variable box length L. With regard to the electron theo- 
retical description of the interatomic potentials in tran- 
sition metal alloys by Hausleitner and Hafner E3' we 
model the interatomic couplings as in [261 ] by a volume 
dependent electron-gas term E vo i (V) and pair potentials 
4>{r) adapted to the equilibrium distance, depth, width, 
and zero of the Hausleitner-Hafner potentials [3Cj for 
Ni .2Zr .8 E3- For this model simulations were started 
through heating a starting configuration up to 2000 K 
which leads to a homogeneous liquid state. The system 
then is cooled continuously to various annealing temper- 
atures with cooling rate —d{T = 1.5xl0 12 K/s. After- 
wards the obtained configurations at various annealing 
temperatures (here 1500-800 K) are relaxed by carrying 
out additional isothermal annealing run at the selected 
temperature. Finally the time evolution of these relaxed 
configurations is modelled and analyzed. More details of 
the simulations are given in |3l| . 



The normalized partial- and incoherent intermediate 
scattering functions are given by 



= F ij {q ) t)/S ij {q) 
$i(q,t) = F*(q,t) , 



(4) 
(5) 



where the Sij(q) — Fij(q,t — 0) are the partial static 
structure factors. 

The basic equations of the MCT are the set of nonlin- 
ear matrix integrodifferential equations given by 



t) = 



F(q,t) + n 2 (q)F(q,t)+ f drM(q,t - r)F(q, 
Jo 

(6) 

where F is the 2x2 matrix consisting of the partial inter- 
mediate scattering functions F i j{q 1 t), and the frequency 
matrix Jl is given by 



[tf(q)]„ 



q' 



(7) 



S(q) denotes the 2x2 matrix of the partial structure 
factors Sij(q), Xi = Ni/N and m; means the atomic mass 
of the species i. 

The MCT for the idealized glass transition predicts @ 
that the memory kern M can be expressed at long times 
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by 



My(q,i) 



k B T 



dk 



EE 



2pm i x j J (2tt) h kiy 

xV i ki(q.,k)Vjk'i'(q., q- k) 
xF fcfc /(k,t)fi,/(q-k,t) , (8) 

where p = N/V is the particle density and the vertex 
Via/3 (q,k) is given by 

^(q,k) = ^(W(k) + q,(q ' k) faca (q - k) (9) 

q q 

and the matrix of the direct correlation function is de- 
fined by 



cy(q) 



[s- x (q) 



(10) 



The equation of motion for F 1 / (g, i) has a similar form 
as eq. ®, but the memory function for the incoherent 
intermediate scattering function is given by: 



M/(q,f) 



dk 1 /qk 



(2tt) 3 p V 9 
xi?(q-k,t), 



(cF)<(M) 



(11) 



(cF)i(k,t) = {c ii {q)) 2 F ii {q,t) + 2c lt {q)c t] {q)F ij (q,t) 



+ (c i j(q)) 2 F jj (q,t) j^i 



(12) 



In order to characterize the long time behaviour of the 
intermediate scattering function, the non-ergodicity pa- 
rameters f(q) are introduced as 



fij(q) = lim t ^ 00 ^ ij (q > t) 



(13) 



These parameters are the solution of eqs. I^J-JTUJ at long 
times. The meaning of these parameters is the following: 
if /ij(q) = 0, then the system is in a liquid state with 
density fluctuation correlations decaying at long times. 
If fij(q) > 0, the system is in an arrested, nonergodic 
state, where density fluctuation correlations are stable 
for all times. In order to compute /r/(q), one can use the 
following iterative procedure |22| : 



f(l+D( a] - A (g)+ B (g) 

W C(q)+D(q) ' 



(14) 



where the matrix A(q), T5(q),C(q), T)(q), N(g) is given 
by 



A(g) = S(g).N[f»f«]( ? ).S(«/) 



(15) 



D(g) = g- 2 |S(g)||N[f«,f«](g)| 



x t k B T 



M t3 (q) 



(18) 



(19) 



This iterative procedure, indeed, has two type of solu- 
tions, nontrivial ones with f(q) > and trivial solutions 
%) = 0. 

The incoherent non-ergodicity parameter //(g) can be 
evaluated by the following iterative procedure: 



: Z"' +1 (g) 
l-ft l+1 (l) 



M?[ZJ S A{q) 



(20) 



As indicated by ea. (|20|l . computation of the incoher- 
ent non-ergodicity parameter f" (q) demands that the co- 
herent non-ergodicity parameters are determined in ad- 
vance. 



IV. RESULTS AND DISCUSSIONS 

A. Partial structure factors and intermediate 
scattering functions 

First we show the results of our simulations concerning 
the static properties of the system in terms of the partial 
structure factors Sij (q) and partial correlation functions 

To compute the partial structure factors S ii(q ) for a 
binary system we use the following definition 35] 

Sij(~cf ) = %&j + P x i x j f(9ij(r) - l)e~ 4 q ' r d~r> {2}) 
where 




are the partial pair correlation functions. 

The MD simulations yield a periodic repetition of the 
atomic distributions with periodicity length L. Trunca- 
tion of the Fourier integral in Eq. I|21(l leads to an oscilla- 
tory behavior of the partial structure factors at small q. 
In order to reduce the effects of this truncation, we com- 
pute from Eq. (|22|l the partial pair correlation functions 
for distance r up to R c = 3/2L. For numerical evaluation 
of ea. (|21J) . a Gaussian type damping term is included 



B( 9 ) =(Z - 2 |S( g )||N[f«,fW]( 9 )|S( 9 )| , (16) 
C( q )=q 2 +Tr(S(q)-*N[{V\f^](q)) , (17) 



Sij(q) = XiSij + AnpxiXj I r 2 (g lj (r) - 1) 



sin(qr) 



qr 



x exp(-(r/i?) 2 )dr 



(23) 
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Q(A" 1 ) 

FIG. 1: Partial structure factors at T = 1400 K, 1300 K, 
1200 K, 1100 K, 1000 K, 900 K and 800 K (from top to bot- 
tom); a) Ni-Ni-part, the curves are vertically shifted by 0.05 
relative to each other; b) Ni-Zr-part, the curves are vertically 
shifted by 0.1 relative to each other; and c) Zr-Zr-part, the 
curves are vertically shifted by 0.5 relative to each other. 

with R = R c /3. 

FigE] shows the partial structure factors Sij(q) versus 
q for all temperatures investigated. The figure indicates 
that the shape of Sij (q) depends weakly on temperature 
only and that, in particular, the positions of the first 
maximum and the first minimum in Sij(q) are more or 
less temperature independent. 

In order to compare our calculated structure factors 
with experimental ones, we have determined the Faber- 
Ziman partial structure factors dij(q) |37j ] 

= l + pJ fa - ly-^-^tr , (24) 

and the Faber-Ziman total structure factor Sf a f(q) [3t| . 
For a binary system with coherent scattering length bi of 
species i the following relationship holds: 

S tJ(9) = 7772 [ x l b la u (q) + x\b 2 2 a 22 (q) 
W 

+2xix 2 bib 2 ai 2 (q)] . (25) 

In the evaluation of a« (q) , we applied the same algo- 
rithm as for Sij(q). By using cnj(q) and with aids of the 
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FIG. 2: Comparison between our MD-simulations and exper- 
imental results [%| of the total Faber-Ziman structure factor 
Stot (<?) and the partial Faber-Ziman structur factors atj(q) 
for Nio.2Zro.8- 

experimental data of the average scattering length b one 
can compute the total structure factor. Here we take bi 
from the experimental data of Kuschke ■ b for natural 
Ni is 1.03 (1(T 12 cm) and for Zr 0.716 (lO^ 12 cm). FigEl 
compares the results of our simulations with the experi- 
mental results by Kuschke [23( for the same alloy system 
at 1000 K. There is a good agreement between the exper- 
imental and the simulations results which demonstrates 
that our model is able to reproduce the steric relations 
of the considered system and the chemical order, as far 
is visible in the partial structure factors. 



B. Non-ergodicity parameters 

The non-ergodicity parameters are defined over 
Eci. (|13|) as a non- vanishing asymptotic solution of the 
MCT-eq. ||BJ • Phenomenologically, they can be estimated 
by creating a master curve from the intermediate scat- 
tering functions with fixed scattering vector q at different 
temperatures. The master curves are obtained by plot- 
ting the scattering functions Q(q, t) as function of the 
normalized time t/r a . 

Fig. |31 presents the estimated g-dependent non- 
ergodicity parameters from the coherent scattering func- 
tions of Ni and Zr, Fig.0]those from the incoherent scat- 
tering functions. In Fig.0aiid0]are also included the de- 
duced Kohlrausch- Williams- Watts amplitudes A(q) from 
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the master curves and from the intermediate scattering 
functions at T=1100 K. (The further fit-parameters can 
be found in |3lj.) 

In order to compute the non-ergodicity parameters 
fij(q) analytically, we followed for our binary system 
the self-consistent method as formulated by Nauroth and 
Kob an d as sketched in Section III. A. Input data for 
our iterative determination of fij(q) = Fij(q, 00) are the 
temperature dependent partial structure factors Sij(q) 
from the previous subsection. The iteration is started 
by arbitrarily setting PVt-jv* (<?, 00) W = 0.55jvi-jvi(9)) 
Fjr r _ Zr (g,oo)W = 0.55 Zr _ Zr (g), F Ni - Zr (q,oo)W = 0. 




q (A" 1 ) 

FIG. 3: Non-ergodicity parameter fdj for the coherent in- 
termediate scattering functions as solutions of eqs. (JJ and 
©(solid line), KWW-parameter A(q) of the master curves 
(diamond), Von Schweidler-parameter f c (q) of the master 
curves (square), and KWW-parameter A(q) for $ij(q) at 1100 
K (triangle up); a) Ni-Ni-part and b) Zr-Zr-part. 

For T > 1200 K we always obtain the trivial solu- 
tion fij(q) = while at T = 1100 K and below we get 
stable non- vanishing fij(q) > 0. The stability of the 
non-vanishing solutions was tested for more than 3000 
iteration steps. From this results we expect that T c for 
our system lies between 1100 and 1200 K. To estimate 
T c more precisely, we interpolated Sij(q) from our MD 
data for temperatures between 1100 and 1200 K by use 
of the algorithm of Press et.al. 39]. We observe that 
at T = 1102 K a non-trivial solution of fij(q) can be 
found, but not at T = 1105 K and above. It means that 
the critical temperature T c for our system is around 1102 
K. The non-trivial solutions fij(q) for this temperature 
shall be denoted the critical non-ergodicity parameters 
fdj(q)- They are included in Fig. [3] As can be seen 
from Fig. |3 the absolute values and the g-dependence 




q (A" 1 ) 

FIG. 4: The same as fig^lbut for the incoherent intermediate 
scattering function; a) Ni-part and b) Zr-part. 

of the calculated f c ij(q) agree rather well with the esti- 
mates from the scattering functions master curve and, in 
particular, with the deduced Kohlrausch- Williams- Watts 
amplitudes A{q) at 1100 K. 

By use of the critical non-ergodicity parameters fdj (q) , 
the computational procedure was run to determine the 
critical non-ergodicity parameters /^(g) for the incoher- 
ent scattering functions at T = 1102 K . Fig. 0] presents 
our results for so calculated /^(g). Like Fig. |3] for the 
coherent non-ergodicity parameters, Fig. ^demonstrates 
for the fci(q) that they agree well with the estimates from 
the incoherent scattering functions master curve and, in 
particular, with the deduced Kohlrausch- Williams- Watts 
amplitudes A(q) at 1100 K. 

C. Diffusion-coeffient 

From the simulated atomic motions in the computer 
experiments, the diffusion coefficients of the Ni and Zr 
species can be determined as the slope of the atomic mean 
square displacements in the asymptotic long-time limit 

(1/NA £ K(t) - r a (0)\ 2 
A(T)=lim . (26) 

t^oo tot 

Fig-Elshows the thus calculated diffusion coefficients of 
our Nio.2Zro.8 model for the temperature range between 
800 and 2000 K. At temperatures above approximately 
1250 K, the diffusion coefficients for both species run par- 
allel with temperature in the Arrhenius plot, indicating 
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FIG. 5: Diffusion coefficients Di as a function of 1000/T. 
Symbols are MD results for Ni (square) and Zr (diamond); 
the full line are a power-law approximation for Ni and for Zr. 
resp.. 

a fixed ratio Dpjil ' Dz r ~ 2.5 in this temperature regime. 
At lower temperatures, the Zr atoms have a lower mobil- 
ity than the Ni atoms, yielding around 900 K a value of 
about 10 for D^i/Dzr- That means, here the Ni atoms 



carry out a rather rapid motion within a relative immo- 
bile Zr matrix. 

According to the MCT, above T c the diffusion coeffi- 
cients follow a critical power law 

Di(T) ~ (T - T c ) 7 , for T > T c (27) 

with non-universal exponent 7 0, l38|| . In order to es- 
timate T c from this relationship, we have adapted the 
critical power law by a least mean squares fit to the sim- 
ulated diffusion data for 1050 K and above. The results 
of the fit are included in Fig. by dashed lines. Ac- 
cording to this fit, the system has a critical temperature 
of 950 K. The parameters 7 turn out as 1.8 for the Ni 
subsystem and 2.0 for the Zr system. 

V. CONCLUSION 

The results of our MD-simulations show that our sys- 
tem behave so as predicted by MCT in the sense that the 
diffusion coefficients follow the critical power law. After 
analizing this coefficient we found that the system has 
critical temperature of 950 K. 

Our analysis of the ergodic region ( T > T c ) and 
of the non-ergodic region ( T < T c ) lead to T c - 
estimations which agree each other within 10 %. These 
T c -estimations are also in acceptable compliance with the 
T c -estimation from the dynamic phenomenons. Within 
the scope of the precision of our analysis, the critical 
temperatur T c of our system is about 1000 K. 
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